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Abstract 

. 1^ The cycle-by-cycle variations in heat release for a simulated spark-ignited engine are 

^ analyzed within a turbulent combustion model in terms of some basic parameters: 

the characteristic length of the unburned eddies entrained within the flame front, a 
characteristic turbulent speed, and the location of the ignition kernel. The evolution 
of the simulated time series with the fuel-air equivalence ratio, (p, from lean mixtures 
to over stoichiometric conditions, is examined and compared with previous experi- 
ments. Fluctuations on the characteristic length of unburned eddies are found to be 
essential to simulate heat release cycle-to-cycle variations and recover experimental 
results. Relative to the non-linear analysis of the system, it is remarkable that at 
fuel ratios around ~ 0.65, embedding and surrogate procedures show that the 
dimensionality of the system is small. 
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1 Introduction 



Nowadays, primary objectives for vehicle manufacturers are the reduction of 
fuel consumption as well as the reduction of emissions. Among the several 
solutions to be adopted are fast combustion, lean burn, variable valve timing, 
gasoline direct injection and some others fi\. Combustion in spark-ignition 
engines in lean burn conditions can lead to cycle-to-cycle variations (CV) 
that are responsible for higher emissions and limit the practical levels of lean- 
fueling. Better models for simulating the dynamic character of CV could lead 
to new design and control alternatives for lowering emissions and enlarging 
the actual operation limits. 

The cyclic variation experienced by spark ignition engines is essentially a com- 
bustion problem which is affected by many engine and operating variables like 
fuel properties, mixture composition near spark plug, charge homogeneity, ig- 
nition, exhaust dilution, etc [2] . Physical sources of CV arise from the merging 
of at least the following elements [31H] : 

(1) Gases mixture motion and turbulences during combustion. 

(2) The amounts of fuel, air, and residual or recirculated exhaust gas in the 
cylinder. 

(3) Homogeneity of the mixture composition, specially in the vicinity of the 
spark plug. 

(4) Details of the spark discharge, like breakdown energy and the initial flame 
kernel random position. 

Not all these factors are equally important, and some others were not included 
in the list. It is usually argued that the first item in the list is the leading effect. 
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The flow motion in any piston engine is unsteady, having an average pattern 
and a fluctuating velocity component. Nevertheless, it is not easy to select an 
adequate turbulence parameter. During the engine operation not all the burnt 
gases are expelled from the combustion chamber in the exhaust stroke. Some 
fraction is mixed with the next incoming air and fuel mixture in the following 
intake stroke. So, when combustion approaches the lean flammability limit, 
very small modiflcations due to recycled gases in the mixture composition or 
temperature can greatly affect the next combustion event causing a highly 
non-linear cycle-to-cycle coupling process. 

Besides experimental and thermodynamic theoretical studies [5], computer 
simulations [6l|7] are useful tools which provide complementary understanding 
of the complex physical mechanisms involved in the operation of a real engine. 
Several previous studies PH] tried to reproduce experimental CV by incorpo- 
rating explicit models for ignition, combustion chemistry, flame evolution, and 
turbulence. In these models some input parameters are varied in each cycle in 
such a way that cycle-to-cycle combustion also changes, and so the pressure 
and temperature in the cylinder and subsequently heat release, power output 
and efficiency. Abdi Aghdam et al. [4j reproduced satisfactory fits to pressure 
and flame radius data by performing quasi- dimensional computer simulations 
and acting on the rms turbulent velocity. 

In regard to CV phenomenology we mention the following experimental and 
theoretical works, many of them related to the non-linear dynamics involved 
in these variations. Letelier et al. [9] reported experimental results for the vari- 
ation of the in-cylinder pressure versus crank angle for a four-cylinder spark 
heat engine. By reconstructing the phase space and building the Poincare sec- 
tions it was concluded that CV is not governed by a chaotic process but rather 
by a superimposition of a non-linear deterministic dynamics with a stochastic 
component. In the works by Daw et al. [TOlITT] . a physically oriented simple 
mathematical model was reported for the spark ignited engine. The fuel/air 
mass ratio in a cycle was modeled as stochastic through random noise and non- 
linear deterministic coupling between consecutive engine cycles. The predicted 
CV trends for heat release versus fuel-air ratio were compared with experi- 
mental results by analyzing bifurcation plots and return maps for different 
parameters of the model (the mean fuel-air ratio, its standard deviation, and 
the fraction of the un-reacted air and fuel remaining in the cylinder). Also, an 
analysis of the temporal irreversibility in time-series measurements was pre- 
sented by using a symbolic approach to characterize the noisy dynamics 



Litak and coworkers [T3|14|ll5j have reported extensive work relative to cycle- 
to-cycle variations of peak pressure, peak pressure location and heat release 
in a four-cylinder spark ignition engine in terms on spark advance angle and 
different torque loadings under near stoichiometric conditions. The observed 
qualitative change in combustion was analyzed using different statistical meth- 



ods as return maps, recurrence plots [T3|T6] , correlation coarse-grained entropy 
\t7\ and, recently, multi-fractal techniques P^HTH] . and the so-called multi-scale 
entropy (MSE) (or sample entropy) [12]: an improved statistical tool to ac- 
count for complexity measure in time series with multiple temporal or spatial 
scales. 

Along this line, a main goal of this paper is to analyze the effect of some com- 
bustion parameters and their consequences on cycle-to-cycle variations of heat 
release for fuel-air equivalence ratios ranging from very lean mixtures to stoi- 
chiometric conditions. To get this we developed a quasi-dimensional computer 
simulation that incorporates turbulent combustion and a detailed analysis of 
the involved chemical reaction that includes in a natural way the presence of 
recycled exhaust gases in the air-fuel fresh mixture and valves overlapping. 
The model, previously developed and validated j^, allows for a systematic 
study of the influence on CV of three basic combustion parameters: the char- 
acteristic length of the unburned eddies, the characteristic turbulent speed, 
and the location of the ignition site. Our simulations correctly reproduce the 
basic statistical parameters for the fuel ratio dependence of experimental heat 
release time series [TTfTH] . Moreover, we shall employ some usual techniques 
from non-linear data analysis in order to characterize our results and compare 
them with those from other authors: first-return maps, bifurcation-like plots, 
and embedding and surrogate techniques to obtain the correlation dimension. 

The paper is organized in the following way. Section II is devoted to briefly 
describe the basic elements of our quasi-dimensional computer simulation with 
special emphasis on combustion. In Sec. Ill we analyze the simulation results 
for heat release time series at different values of the fuel/air equivalence ratio 
and compare with previous experimental results. Sec. IV is devoted to ana- 
lyze the dimensionality of the system by means of embedding and surrogate 
procedures. Finally, in Sec. V we discuss our results and summarize the main 
conclusions of the paper. 



2 Simulated model: basic framework 



We make use of a quasi- dimensional numerical simulation of a monocylin- 
drical spark ignition engine, previously described and validated [6]. Opposite 
to zero-dimensional (or thermodynamical) models where all thermodynamic 
variables are averaged over a finite volume and an empirical correlation is used 
to approximate the combustion process, within the quasi-dimensional scheme 
explicit differential equations are considered to solve combustion under the 
hypothesis of a spherical flame front. In our simulation model two ordinary 
differential equations are solved in each time step (or each crankshaft angle) 
for the pressure and the temperature of the gases inside the cylinder. This set 



of coupled equations (that are explicitly written in, [6] , [7], and [20]) are 
valid for any of the different steps of the engine evolution (intake, compres- 
sion, power stroke, and exhaust) except combustion. For this stage a two- zone 
model discerning between unburned [u) and burned (6) gases that are sepa- 
rated by an adiabatic flame front with negligible volume is considered. Thus, 
the temperature equation splits in equations for T^ and T^. Moreover, differen- 
tial equations giving the evolution of the unburned and burned mass of gases 
should be solved. Next we detail the main assumptions of the combustion 
model considered. 



2. 1 Combustion process 



Combustion is the most important process in cycle-by-cycle variations. For 
simulating combustion we assume the turbulent quasi-dimensional model (some- 
times called eddy- burning or entrainment model) proposed by Blizard and 
Keck prp^ and improved by Beretta et al. [23j, by considering a diffusion 
term in the burned mass equation as we shall see next. The model starts from 
the idea that during flame propagation not all the mass inside the approx- 
imately spherical flame front is burned, but there exist unburned eddies of 
typical length It. Thus, a set of coupled differential equations for the time evo- 
lution of total mass within the flame front, rrie (unburned eddies plus burned 
gas) and burned mass, mb, is solved: 

rrie-mb , . 

rrib = PuAfSi H (1) 

The = PuAf [ut (l - e*/"^) + Si\ (2) 

where Ut is a characteristic velocity at which unburned gases pass through the 
flame front, pu is the unburned gas density, and Aj is the flame front area, r^ 
is a characteristic time for the combustion of the entrained eddies, r;, = U/ Si 
and Si is the laminar combustion speed that is determined from its reference 
value, SiQ, as [24j: 




where Hr is the mole fraction of residual gases in the chamber and a and /3 are 
functions of the fuel/air equivalence ratio, (j) [2S]. The reference laminar burn- 
ing speed. Sift, at reference conditions (Tref, Pref) is obtained from Giilder's 
model [26|25] . Because of the reasons that we shall detail later, it is important 
to stress that the laminar speed, apart from the thermodynamic conditions, 
depends on the fuel/air equivalence ratio and on the mole fraction of gases 
in the chamber after combustion. Thus, Eq. ^ shows the coupling between 



combustion dynamics and the residual gases in the cyhnder after the previ- 
ous combustion event. In other words laminar burning speed depends on the 
memory of the chemistry of combustion. From Eq. ^ it is easy to see that 
Tb determines the time scale associated either to the laminar diffusive evolu- 
tion of the combustion at a velocity Si or to the rapid convective turbulent 
component associated to Ut- 

In order to numerically solve the differential equations for masses (Eqs. ([I]) 
and (|2])), it is necessary to determine Aj, Ut, and It- The flame front area, 
Af, is calculated from the flame front radius, considering a spherical flame 
propagation in a disc shaped combustion chamber [22]. The radius is calcu- 
lated through the enflamed volume, Vf, following the procedure described by 
Bayraktar |27] . 

TTlb , TTie-mb ,... 

Vf = \ [4:) 

Pb Pu 

The flame front area depends on the relative position of the kernel of com- 
bustion respect to the cylinder center. Re, which have its nominal value at 
the spark plug position, but convection at early times can produce signifi- 
cant displacements. For Ut and It, Beretta et al. [23] have derived empirical 
correlations in terms of the ratio between the fresh mixture density at ambi- 
ent conditions, pi, and the density of the unburned gases mixture inside the 
combustion chamber, p„: 




Mi = 0.08MJ^ (5) 



/i = 0.8L,,^axf-j' (6) 

where Ui is the mean inlet gas speed, and L^^max the maximum valve lift. 
Equations ([I]) and ^ must be solved simultaneously together with the set of 
differential equations for temperature and pressure [6]. Respect to the end of 
combustion, we assume that if exhaust valve opens before the completion of 
combustion, it finishes at the moment of the opening. 

Among the parameters that determine the development of combustion, there 
are three essential ones: the characteristic length of eddies U (associated to 
the characteristic time, r^), the turbulent entrainment velocity ut (essential 
for the slope of me{t) during the fastest stage of combustion), and also Af or 
the flame radius center Rf. that gives the size and geometry of the flame front 
(and thus influences all process). Therefore, all these parameters could have 
strong influence on cycle-by-cycle variations. 

In order to calculate the heat release during combustion we apply the first 
principle of thermodynamics for open systems separating heat release, 5Qr, 
internal energy variations associated to temperature changes, dU , work out- 



put, SW, and heat transfers from the working fluid (considered as a mixture 
of ideal gases) to the cylinder walls, 6Qe: 

5Qr = dU + 6W + 5Qi (7) 

where internal energy and heat losses include terms associated to either un- 
burned or burned gases: U = muC^^uTu + 'mbCvfiTb and SQ^ = 6Qi^u + SQi,b- 
All these terms can be derived in terms of time or of the crankshaft angle. 
Net heat release during the whole combustion period is calculated from the 
integration of heat release variation during that period. A particular model to 
evaluate heat transfer through cylinder walls should be chosen before the nu- 
merical obtention of Qr- Units throughout the paper are international system 
of units except when explicitly mentioned. 



2.2 Heat transfer 



Among the different empirical models that can be found in the literature for 
heat transfer between the gas and cylinder internal wall we shall utilize that 
developed by Woschni [28] . 

Qi 1 or, o „,0.8„,,0.8 D-0.2Tn-0.55 



At (T - T^) 



129.8 p^-^'w^-^B-^-'T-"-'' (8) 



where At is the instantaneous heat transfer area, p is the pressure inside the 
cylinder, T the gas temperature, T^u the cylinder internal wall temperature, B 
the cylinder bore, and w the corrected mean piston speed, that is calculated 
as follows, 

W = CiVp + C2^ri^ip-Pm) (9) 

PrVr 

Ci and C2 are constant parameters |i29j, vp is the mean piston speed, Vdt is the 
maximum displaced volume, pm is a motoring pressure [29], and the subscript 
r refers to conditions immediately after inlet valve is closed. All units are in 
the international system of units, except pressures which are in bar. 



2.3 Chemical reaction 



In order to solve the chemistry and the energetics of combustion we assume 
the unburned gas mixture as formed by a standard fuel for spark ignition 
engines (iso-octane, CgHig), air, and exhaust gases. Exhaust chemical compo- 
sition appears directly as calculated from solving combustion. The considered 



chemical reaction is, 

(1 - Ur) [CgHis + a(02 + 3.773N2)] + 

Vr /3rC02 + 7rH20 + /i^N2 + Z/^02 + SrCO + (5^2 + 

+ ^^H + K^O + cr^OH + ^30 

/3CO2 + 7H2O + /iN2 + z/02 + eCO + ^2 



12- 



+ ^H + kO + aOH + gNO (10) 

where the subscript r refers to residual. We make use of the subroutine de- 
veloped by Ferguson [3Dj (but including residual gases among the reactants) 
to solve combustion and calculate exhaust composition. Our model does not 
consider traces of CsHig in combustion products, but in the energy release we 
actually take into account combustible elements as CO, H, or H2. The thermo- 
dynamic properties of all the involved chemical species are obtained from the 
constant pressure specific heats, that are taken as 7-parameter temperature 
polynomials [3Tj . 



As a summary of the theoretical section it is important to note that in our 
dynamical system the coupled ordinary differential equations for pressure and 
temperature are in turn coupled with two other ordinary differential equations 
for the evolution of the masses during combustion. So, globally this leads to a 
system of differential equations where apart from several parameters (mainly 
arising from the geometry of the cylinder, the chemical reaction and the models 
considered for other processes), there is a large number of variables: pressure 
inside the cylinder, p, temperatures of the unburned and burned gases, T^, 
and Tfc, and masses of the unburned and burned gases, niu and m^. All these 
variables evolve with time or with the crankshaft angle. So, up to this point, 
our dynamical model is a quite intricate deterministic system with those time 
dependent variables. 



3 Results 



In this section we present the numerical results obtained from the simulations 
outlined above and taking the numerical values for the cylinder geometry from 
Beretta et al. |23] and for a fixed engine speed of 109 rad/s. Details on some 
running parameters of the computations can be found in [B] . 

The computed results for heat release, Qr, for each of the first 200 cycles 
are presented in Fig. 1(a) for a particular value of the fuel-air equivalence 
ratio, = 1.0. This time evolution was obtained directly from the solution of 



the deterministic set of differential equations. The curve clearly shows, after 
a transitory period, a regular evolution that does not match with previous 
experimental studies of cycle-to-cycle variability on heat release [HHn]- For 
other fuel ratios the curves obtained present different levels of variability but 
never display the typical experimental fluctuations. In order to analyze these 
time series we show in Table I for several fuel ratio values, some usual statistical 
parameters: the average value n, the standard deviation a, the coefficient of 
variance, COV = a/ ft, the skewness S = J^iLii^i ~ t^Y/li^ ~ l)*^^]; and 
the kurtosis K = J2f=iixi — ^)^/[{N — l)cr^]. 5" is a measure of the lack of 
symmetry in such a way that negative (positive) values imply the existence of 
left (right) asymmetric tails longer than the right (left) tail. K is a measure 
of whether the data are peaked {K > 3) or flat {K < 3) relative to a normal 
distribution. 

With the objective to recover the experimental results on cyclic variability we 
have checked the influence of incorporating a stochastic component in any of 
the physically relevant parameters of the combustion model. We first analyze 
the infiuence of the characteristic length It and velocity Ut keeping the location 
of ignition at the spark plug position {i.e., Re = 25 x 10~^ m [23])- Although 
both parameters can be considered as independent in our model (see below), 
we assume that they are linked by the empirical relations (IS]) and (|6]). We 
fitted the experimental results by Beretta [23] for It considered as a random 
variable to a log- normal probability distribution, LogN{fiiogi^,aiogiJ, around 
the nominal value /° = O.SLy^rnaxipi/ pu)^^^ [Eq- (p|] with standard deviation 
ciogit = 0.222 and mean /^logit = log(/°) — {aiggij2). Then, in each cycle Ut is 
obtained from Eq. ([s]) through the empirical densities ratio. 

Table II for several values of ranging from very lean mixtures to over stoi- 
chiometric, show the characteristics of the time evolution of the heat release 
when the stochastic component on It has been introduced. Figure 1(b) displays 
a time series for the stochastic computations at least qualitatively much simi- 
lar to the experimental ones than the deterministic for = 1.0. It is important 
to note that it is difficult to perform a direct quantitative comparison with ex- 
periments, because of the large number of geometric and working parameters 
of the engine [6] , usually not specified to the full extent in the experimental 
publications. 

From a comparison between Tables I and II is clear that average values are 
of course similar, but the standard deviation and covariance are much smaller 
in the deterministic case. The behavior of skewness is subtle, it only could be 
concluded that globally in the deterministic case values for S are closer to zero 
(that would be the value corresponding to a Gaussian distribution). For the 
stochastic simulation kurtosis in the fuel ratio interval between = 0.6 and 
0.8 is much higher than in the deterministic case, that is a sign of more peaked 
distributions. Figure 2 shows the evolution of heat release with the number of 



cycles for fuel ratios between 0.5 and 1.1. A careful inspection shows that at 
low and intermediate fuel-ratios distributions are quite asymmetric with tails 
displaced to the left. In other words there exist remarkable poor combustion 
or misfire events. This is quantified by the evolution with (p of the skewness 
compiled in Table II: it takes high negative values in the interval (f) = 0.5 — 0.8. 
These results are in accordance with the experimental ones by Sen et al. [18] 
(see Fig. 1 in that paper). 

We represent in Fig. 3 the evolution with the fuel-air equivalence ratio of the 
statistical parameters contained in Tables I and II in order to perform a direct 
comparison with previous experimental results. Although the experiments by 
Daw et al. [TTfTS] were performed for a real V8 gasoline engine with a different 
cylinder geometry that the considered in our simulations, the dependence on 
of the standard deviation, covariance, skewness and kurtosis is very similar, 
although, of course, vertical scales in simulations and experiments are differ- 
ent. It makes no sense, to compare the mean values of heat release because 
the different shape and size of the real and the simulated engine. It is clear 
from the figure that at intermediate fuel ratios, around, = 0.7, skewness also 
reproduces a pronounced minimum and kurtosis a sharp maximum. Neverthe- 
less, our deterministic heat release computations do not reproduce neither the 
evolution with shown by experiments nor the magnitude of those statistical 
parameters. 

Figure 4 represents the first-return maps, Qr,i+i versus Qr,ii for the same val- 
ues of (j) obtained from the deterministic simulation (dark points) and also 
from the stochastic computations (color points). Return maps when It does 
not have a stochastic component, for any value of fuel ratio, seem noisy un- 
structured kernels. When It is considered as stochastic a rich variety of shapes 
for return maps is found, depending on the fuel ratio. From the stochastic 
return maps in this figure and the statistical results in Table II we stress the 
following points. First, at high values of (0.9 — 1.1) variations of heat release 
behave as small noisy spots characteristics of small amplitude distributions 
with asymmetric right tails and slightly more peaked near the mean than a 
Gaussian distribution. For = 1.0 and 1.1 those noisy spots are partially 
overlapped. These results are very similar to the experimental ones by Daw 
et al. (see Fig. 1 in |32]) for different engines and those by Litak et al. (see 
Fig. 4 in [I3]). Second, at intermediate fuel ratios (0 = 0.6 — 0.8) extended 
boomerang-shaped patterns are clearly visible. Note that these series lead to 
distributions very peaked near the mean, decline rather rapidly, and have quite 
asymmetric left tails. In this interval our results also are in accordance with 
those of Daw et al. [TT|32] . Third, at the lowest equivalence ratio, = 0.5 a 
change of structure is observed, probably because of cycle misfires for this very 
poor mixture. Here the probability distribution is less peaked than a Gaus- 
sian, present a high standard deviation, and an asymmetric left tail. Fourth, a 
closer inspection of Fig. 4 reveals some kind of asymmetry around the diago- 
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nal, more pronounced as the fuel-air ratio becomes lower. This is in accordance 
both with experimental and model results by Daw et al. [TIHllf32j . So, our 
stochastic scheme with the consideration of pseudo-random fluctuations on It 
reproduces the characteristic return maps of this kind of systems from poor 
mixtures to over stoichiometric ones. 

The results above refer to the the joint influence of the parameters Ut and 
It on the CV-phenomena when both parameters are considered linked by the 
empirical relations Eqs. (pi) and (pi). From now on our goal is double: on one 
side, to analyze the influence of Ut and It but considered as independent in 
the simulation model and, on the other side, to analyze the influence of the 
third key parameter in the combustion description, the location of ignition 
accounted for the parameter R^.. To get this we have run the simulation for 
2800 consecutive cycles for the previously considered fuel-air equivalence ratios 
but taking these three parameters, one by one, as stochastic in nature. 

First, we have considered only fluctuations on It, taking the same distribution 
that at the beginning of this section. We show the corresponding return maps 
in Fig. 5(a), that should be compared with Fig. 4. From this comparison we 
can conclude that the main characteristics of the return maps obtained when 
It and Ut were considered as linked (noisy spots near stoichiometric conditions, 
boomerang-shaped structures at intermediate fuel-air ratios and complex pat- 
terns at lower air-fuel ratios) are already induced by the distribution of the 
characteristic length It when the other parameters are not stochastic. Differ- 
ences between both figures only affect the dispersion of the boomerangs arms. 

In regard to Ut-, it is remarkable that it is not easy to find in the litera- 
ture experimental data which allow to deduce stochastic distributions for this 
characteristic turbulent velocity with certainty. We have used the data in|l] 
for the turbulence intensity and translated them to generate our log-normal 
distribution for the characteristic velocity, Wt, and checked slight changes 
in the distribution parameters within realistic intervals. So, we assume a 
log-normal distribution LogN{n\o^uti^\ogut) around the nominal value m° = 
0.08'Uj(pu/pj)^/^ [Eq. (p|] with standard deviation (J\o^ut = 0.02 m° and mean 
fJ'iogut = log(^t) ~ (^iogMt/2)- When ut is the only stochastic parameter, the 
observed behaviors are not significantly altered (Fig. 5(b)): boomerang-like 
arrangements are found only at low 0, although some sensitivity in the arms 
of the boomerangs is appreciable. Moreover, we have not found in the return 
maps new features respect to those generated with It. 

Respect to the ignition location. Re, some authors [33] concluded from ex- 
periments and models that displacement of the flame kernel during the early 
stages of combustion has a major part in the origination of cycle-by-cycle 
variations in combustion. So, it seems interesting to check this point from 
quasi-dimensional simulations. It is remarkable that it is not straightforward 
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to numerically evaluate the radius and area of the flame front when R^ is 
displaced respect to the origin. We have generalized the studies by Bayrak- 
tar 133] and Blizard [22] to obtain numerical expressions for the radius and 
area of the flame front for whichever value of Re and at any time during flame 
evolution. Second, we assume a Gaussian distribution with standard deviation 
aR^ = 2.965 X 10^'^ m, according to the results reported by Beretta [23]. At 
sight of Fig. 5(c) it is clear that the influence of a stochastic component in Re 
is limited to cause noisy spots at any fuel ratio, clear structures as boomerangs 
are not found at any fuel-air equivalence ratio. So, fluctuations on the displace- 
ment of the initial flame kernel seem not enough to reproduce characteristic 
patterns of CV. 

Finally, we have checked what happens when the three parameters are si- 
multaneously and independently introduced as stochastic in the simulations 
with the distributions mentioned above and no new features were discovered. 
At sight of these results we can say that the observed heat release behavior 
in terms of the fuel-air ratio is mainly due to stochastic and non-additive 
variations of It and Ut and the non-linearity induced by the dynamics of the 
combustion process, independently of ignition location fluctuations. 

Some physical consequences of cyclic variability are susceptible of a clear in- 
terpretation in terms of the evolution of the fraction of burned gases in the 
chamber with time or with crank angle. This evolution has two main ingredi- 
ents: the ignition delay (the interval until this fraction departs from zero) [9] 
and the slope of that increase during combustion (see, for instance. Fig. 16 
in [23]). These factors are mainly associated to eddies characteristic length If 
and Ut, and the laminar speed Si (that in turn is a function of the fraction 
of residual gases, yr, in the chamber, Eq. ([3])). This last variable determines 
the memory effects from cycle to cycle, as noted by Daw et al. [TOfTT] . Cyclic 
variability provokes cycle-to-cycle changes in the shape of the rate of burned 
gases, and this directly affects the evolution of pressure during the cycle, par- 
ticularly its maximum value and its corresponding position. And consequently 
this leads to changes in the performance of the engine, i. e., in the power output 
and in the efficiency. 



4 Nonlinear analysis: Correlation integral and surrogates 



According with the conclusions of the preceding section from now on all the re- 
sults we present were obtained considering fluctuations only on k, except when 
simulations were purely deterministic. In order to get a better comprehension 
of the complexity of heat release fluctuations for different fuel ratio values, 
we next explore the nonlinear properties of the signals through the analysis of 
the correlation integral [5^361137] . The correlation dimension analysis quan- 
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tifies the self-similarity properties of a given sequence and is an important 
statistical tool which helps to evaluate the presence of determinism in the sig- 
nal. We first reconstruct an auxiliary phase space by an embedding procedure. 



X 



J 



The correlation sum C{e,m) = 2/[N{N - 1)] Ei=i E}=i+i Q{e- \\ a 
where B is the Heaviside function, e is a distance and Xk are m-dimensional 
delay vectors, is computed for several values of m and e |38]. If C{e, m) oc e~'^, 
the correlation dimension is defined as d{e) = d\ogC{e,m)/d\oge. To get a 
good estimation of m we need to repeat the calculations for several values of m, 
the number of embedding dimensions. In this way we are able to select which 
one corresponds to the best option to characterize the dimensionality of the 
system. We notice that the application of correlation dimension analysis alone 
is not sufficient to differentiate between determinism and stochastic dynam- 
ics. To reinforce our calculations of correlation integral statistics, we consider 
a surrogate data method to verify that low values of correlation dimension 
of sequences are not a simple consequence of artifacts [39]. To this end, two 
surrogate sets were considered. First, the original sequence was randomly shuf- 
fled to destroy temporal correlations. A second surrogate set was constructed 
by a phase randomization of the original sequence. The phase randomization 
was performed after applying the Fast Fourier Transform (FFT) algorithm to 
the original time series to obtain the amplitudes. Then the surrogate set was 
obtained by applying the inverse FFT procedure [39^0] . 



We apply the correlation dimension method to heat release sequences with 10^ 
cycles and for different values of the fuel ratio in the interval 0.5 < < 0.7. 
For larger fuel ratios we disregard the correlation dimension analysis since the 
signal shows small amplitude fluctuations around a stable value. According 
with the map-like characteristics observed in Fig. 4 we use a delay time equal 
to one. Figure 6 shows the correlation sum for two selected values of fuel 
ratio. We observe that for low values oi (j) {(j) = 0.5, Figs. 6(a)) the correlation 
integral shows a power law behavior against e for several values of m. This 
behavior is best found by plotting the slope d{e) of log(£:, m) vs. loge. The 
inset shows this plot for values of m in the range 1 — 5. In this case there is 
a wide plateau of e which corresponds to the power law behavior. It is also 
clear that the height of the plateau increases with the embedding dimension 
with a reduction of amplitude of the power law region and some fluctuations 
for low values of e. 

In contrast, we observe in Fig. 6(b) that for intermediate fuel ratios (0 = 0.65) 
the scaling region only exists for values of e in the range 10 < e < 10^. This 
is conflrmed by the presence of a plateau for several values of m (see the 
inset in Fig. 5(b)). Remarkably, the height of the plateau, d{e), almost does 
not change with the embedding dimension, indicating that the system can 
be characterized by a very low dimension. We also notice that for very small 
scales {e < 10), a power law behavior can be also identifled, but the scaling 
exponent increases to reach the value of the embedding dimension, suggesting 
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to identify this region as a noisy regime. 

Next, we use the surrogate data method to evaluate if the low value of cor- 
relation dimension observed for cf) = 0.65 is not an artifact. We generated a 
randomly shuffled sequence and another set by a phase randomization of the 
original data. For both surrogate sets the correlation dimension was calcu- 
lated in the same form and intervals as for original data. The results of these 
calculations are presented in Fig. 7. We observe that for shuffled and phase 
randomized data the correlation dimension increases as the value of m also 
increases, indicating that there is no saturation for correlation dimension (Fig. 
7(a) and 7(b)). 

The findings about correlation dimension of original and surrogate data are 
summarized in the plots of Fig. 7(c). For low fuel ratio values, the correla- 
tion dimension does not saturate for high embedding values whereas for the 
intermediate fuel ratio (0 = 0.65), correlation dimension saturates with the 
embedding dimension which is an indication of low dimensional dynamics. It 
is also interesting that when fluctuations on It are not included in the sim- 
ulations (deterministic case) the same situation is observed. The results for 
surrogate data (shuffled and phase randomized) reveal that no saturation of 
correlation dimension is observed, indicating that surrogate sequences differ 
from the original data and that the low dimensionality is probably related to 
the presence of determinism in heat release fluctuations. In this sense, SchoU 
and Russ [H] have reported that deterministic patterns of CV are the conse- 
quence of incomplete combustion, which occurs at = 0.65 in agreement with 
Heywood [2]. 

We present in Fig. 8 an extensive analysis of the evolution of the heat release 
with the fuel ratio between = 0.5 and 1.1. The figure contains the results 
from the direct solution of the deterministic set of equations (black dots) and 
also from incorporating a stochastic component in It (gray dots). For each 
value of 0, the last 100 simulated points of 250 cycles runs are shown. We 
remark that under moderately lean fuel conditions (0 around 0.65) our quasi- 
dimensional stochastic simulation does not show any period-2 bifurcation, as 
it happens in the theoretical model by Daw et al. [TOlITT] . Our work predicts a 
very low dimensionality map in this region. Such conclusion makes sense since 
bifurcations have never been reported in real car engines. 

On the other hand, the deterministic simulations do not lead to a fixed point as 
in the Daw's map, but to a variety of multiperiodic behavior for different fuel 
ratios. In particular, the inset of Fig. 8 shows the evolution of the normalized 
deterministic heat release, Qr/Qn {Qr is the average heat release for each 
value of 0) with fuel ratio between 0.5 and 1.1. We observe that for low and 
high fuel ratio values the number of fixed points is clearly smaller than those 
present for intermediate values of 0, indicating a high multiperiodicity. A 
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deeper evaluation of the multiperiodicity observed in bifurcation-like plots for 
this system probably will deserve future studies. 



5 Discussion and conclusions 



We have developed a simulation scheme in order to reproduce the experi- 
mentally observed fluctuations of heat release in an Otto engine. The model 
relies on the first principle of thermodynamics applied to open systems that 
allows to build up a set of first order differential equations for pressure and 
temperature inside the cylinder. Our model includes a detailed chemistry of 
combustion that incorporates the presence of exhaust gases in the fresh mix- 
ture of the following cycle. The laminar combustion speed depends on the 
mole fraction of residual gases in the chamber through a phenomenological 
law pS]. Combustion requires a particular model giving the evolution of the 
masses of unburned and burned gases in the chamber as a function of time. We 
consider a model where inside the approximately spherical flame front there 
are unburned eddies of typical length It. The model incorporates two other 
parameters, the position of the kernel of combustion, Re, and a convective 
characteristic velocity Ut (related with turbulence), that could be important 
in order to reproduce the observed cycle-to-cycle variations of heat release. 
Actually, we have checked the influence of each of those parameters when 
they are considered as stochastic. It is possible from experimental results in 
the literature to build up log- normal probability distributions for It and Ut, 
and Gaussian for R^. The results of our simulations with these stochastic dis- 
tributions show that the consideration of k or Ut as stochastic is essential to 
reproduce experiments. These parameters can be considered as independent 
or related through an empirical correlation. Both possibilities lead to similar 
results. On the contrary a stochastic component on Re only provokes noisy 
spots in the return maps for any fuel ratio. 

Moreover, an interesting evolution of the temporal series of heat release is ob- 
tained when the fuel ratio, 0, is considered as a parameter. The dependence 
on of the main statistical parameters of heat release sequences for real en- 
gines are reproduced [TTp8] . Specially interesting is the evolution of skewness 
and kurtosis. Return maps also have a rich behavior: from noisy unstructured 
clusters {shotgun patterns) at high fuel ratios to boomerang asymmetric mo- 
tives at intermediate 0. This behavior was obtained before from mathematical 
models and also from experiments [TOlITT] . 

The correlation integral analysis reveals that, for a certain range of the distance 
e, the dimensionality of the system is small only for intermediate fuel ratios. 
We have additionally used the surrogate data method to compare the findings 
of dimensionality between original simulated and randomized data. Our results 
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support the fact that the low dimensionahty observed for intermediate values 
of the fuel-air equivalence ratio (around = 0.65), is related to the presence 
of determinism in heat release fluctuations. When a systematic study of the 
evolution of the heat release calculated from the stochastic simulations with 
is performed, bifurcations were not found. This is in agreement with real 
engines, where to our knowledge period doubling bifurcations have not been 
found. 

In summary, our numerical model incorporates the main physical and chemical 
ingredients necessary to reproduce the evolution with fuel ratio of heat release 
CV. We have particularly analyzed the influence on CV phenomena of three 
basic parameters governing the turbulent combustion process: length of the 
unburned eddies, the turbulent intensity, and the location of the ignition site. 
The influence of the two first ones seems to be basic in the observed structures 
of the heat release time series. The control of these intermittent fluctuations 
of heat release could result in heat engines with greater mean values of power 
and efficiency. Particular rich behavior of heat release at certain low and inter- 
mediate fuel ratios deserves a detailed study of its non-linear dynamics. Work 
along these lines is in progress. 
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Figure 1 : 

Fig. 1. (a) Heat release time evolution as obtained from the simulation without the 
consideration of any stochastic component for (p = 1.0. (b) Heat release time series, 
Qr, obtained with a log- normal distribution (see text) for the characteristic length 
for unburned eddies during combustion, If . 



20 



O) 



1200 ^ 



!W«#?:i^^ 



yw<)l//i((«fMW-^^fM*iil4iMfM^^ 



1000 1 



800 ■ 



600 ■ 



400 ■ 



200 



1.0 



0.9 




0.8 



l^l^t^>f(f«k^^y^fVI(|IH+^'W^^^^ 



Wm 



ffi 





0.6 




1000 



Fig. 2. Heat release time series, Qr, obtained with a log-normal distribution (see 
text) for the characteristic length of unburned eddies during combustion, It- Results 
are shown for several fuel-air equivalence ratio (</> increases from bottom to top from 
6 = 0.5 to 1.1). 
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Figure 3: 

Fig. 3. Evolution of some statistical parameters with the fuel ratio. Left panel, 
experiments |ll|18j and right panel, stochastic simulation results. 
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Figure 4: 

Fig. 4. Heat release return maps after 2800 simulation runs for the same values of (p 
that in Fig. 2 when It is considered as stochastic. For each (p the central dark kernel 
corresponds to the deterministic simulation. 
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Fig. 5. Heat release return maps obtained when the three basic parameters of com- 
bustion, It, ut, and Re are considered as independent and one- by-one stochastic. 
(a) Fluctuations are only introduced in It (see text for details on the probability 
distribution); (b) only ut fluctuates, and (c) only Re is considered as stochastic. 
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Fig. 6. Correlation integral C(e) as a function of the distance e for several embedding 
dimensions m and two values of (j) = 0-5 and (/> = 0.65. The insets show in each case 
logarithmic plots of the correlation dimension vs. e. 
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Figure 8: 

Fig. 8. Evolution of heat release time series as a function of (j). Dark dots correspond 
to the deterministic simulation, and the inhomogeneous gray dots distribution to the 
stochastic approach. The inset shows the normalized heat release, Qr/Qr obtained 
from the deterministic simulations [Q^ is the average heat release for each fuel 
ratio). 
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Table 1 

Statistical parameters of the heat release temporal series, considering the deter- 
ministic model, for several fuel ratio values: mean value, n, standard deviation, a, 
coefficient of covariance, COV, skewness, S, and kurtosis, K. 

0.7 0.8 0.9 1.0 1.1 

886.82 1019.73 1138.79 1227.20 1253.32 

0.785 1.463 1.969 2.249 2.255 

0.885 1.435 1.729 1.832 1.799 

-0.106 -0.013 1.317 0.376 -0.062 

1.332 1.326 3.610 1.361 3.458 

Table 2 

Statistical parameters of the heat release temporal series represented in Fig. 2 and 

obtained considering stochastic fluctuations in If, for several fuel ratio values. 

0.7 0.8 0.9 1.0 1.1 

885.43 1020.24 1139.51 1227.16 1253.71 

17.937 7.333 6.326 6.663 7.093 

20.3 7.2 5.6 5.4 5.7 

-8.388 -3.194 0.147 0.423 0.335 

13.455 104.475 49.491 3.579 3.710 3.454 
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